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Abstract  - In  this  paper  we  present  a simple 
modification  to  the  traditional  Finite  Difference 
Time  Domain  (FDTD)  method  for  treating  material 
cells.  The  Yee  cell  is  split  into  8 smaller  material 
subcells  so  that  each  E and  H field  point  is 
considered  to  be  located  at  the  crosspoint  of  8 cells 
with  differing  material  properties.  Thus  there  is  no 
longer  an  overlap  of  the  material  cells  associated 
with  the  components  of  the  E and  H fields.  The  8 
materia]  properties  are  averaged  at  each  crosspoint. 
Since  the  averaging  is  done  outside  the  time 
marching  loop  there  is  little  increase  in  the  total 
computational  time.  Numerical  results  are  shown  for 
a sinusoidal  plane  wave  scattering  from  a dielectric 
sphere.  These  results  are  compared  with  the  exact 
Mie  solution  and  the  traditional  material  cell  method 
along  different  cuts  through  the  sphere.  The  splitting 
and  averaging  is  shown  to  give  improved  amplitude 
accuracy  in  the  vicinity  of  the  sphere.  Improvement 
is  also  observed  at  planar  interfaces  angled  with 
respect  to  the  grid.  An  additional  benefit  of  this 
subcell  formulation  is  that  objects  may  be  modeled 
with  twice  the  geometrical  resolution  without 
increasing  the  size  of  the  staggered  grid. 

Index  Terms-  FDTD,  Subcell,  Split  cell,  Effective 
Dielectric  Constant 


I.  INTRODUCTION 

The  finite-difference  time  domain  (FDTD) 
technique  is  a widely  used  method  for  computing 
electromagnetic  scattering.  The  formulation  is  based  on 
the  ideas  put  forth  in  [1]  and  involve  computing  the 
electric  and  magnetic  (E  and  H)  fields  on  a staggered 
grid.  There  are  some  shortcomings  with  the  Yee  method 
when  applied  to  interfaces  between  materials.  First,  the 
material  cells  associated  with  the  field  components 
(staggered  in  space)  are  centered  around  each  field  point 
and  thus  overlap  in  space.  Second,  the  discontinuity  in 
permittivity  at  interfaces  between  different  media  causes 
accuracy  problems.  Third,  is  the  error  caused  by 
approximating  a curved  boundary  by  a staircased  one. 

Various  methods  have  been  proposed  to  improve 
the  accuracy  of  the  basic  FDTD  method  at  interfaces, 
mainly  by  reducing  the  staircasing  errors  caused  by 
angled  or  curved  boundaries.  For  example,  one  method 
uses  a volume  weighted  average  of  £ on  each  side  of  a 


material  interface  [2].  An  effective  dielectric  constant  is 
computed  as  follows: 

€eff  = O^i  )/^cctl  > O) 

where  Vx  is  the  volume  on  one  side  of  the  interface  and 

V2  is  the  volume  on  the  other  side  of  the  interface.  This 

method  requires  a special  treatment  of  each  material 
interface.  Specifically  the  method  requires  identifying 
which  cells  are  intersected  by  the  boundary  and  where 
the  boundary  intersects  these  cells.  Dey  and  Mittra 
reported  that  the  results  using  this  approach  are 
significantly  more  accurate  than  using  a staircase  to 
approximate  the  surface  of  a dielectric  structure. 
Another  approach  uses  effective  dielectric  constants 
together  with  the  contour-path  integral  FDTD  (CFDTD) 
method  [3].  Recently  a new  technique  was  developed 
that  utilizes  the  electric  field  components  along  the 
edges  of  the  cell  [4].  Rather  than  using  the  volume  of 
the  dielectric  occupying  a cell  it  instead  uses  the  length 
of  the  dielectric  intersecting  the  cell  edges. 

Still  another  method  uses  dielectric  subgrid 
resolution  (DSR)  [5].  This  can  be  described  as 
partitioning  a standard  cell  into  8 subcells,  each 
homogeneously  filled.  The  8 cells  can  have  different 
material  properties.  Generalized  constitutive  parameters 
are  obtained  by  treating  each  cell  as  an  equivalent 
lumped  circuit.  It  is  shown  by  numerical  example  that 
on  coarser  grids  the  simpler  averaging  method  gives 
more  accuracy  than  the  lumped  circuit  approach  across 
interfaces.  The  DSR  method  also  requires  greater 
memory  since  there  are  effective  dielectric  constants  for 
3 directions. 

The  practical  advantage  to  the  approach  presented 
in  this  paper  is  that  no  special  preprocessing  operations 
are  required  to  compute  the  intersections  of  cells  with 
interface  boundaries.  The  comparisons  in  this  paper  are 
on  the  basis  of  field  amplitudes  as  opposed  to  resonance 
frequencies  used  in  the  cited  articles. 

II.  SIMPLE  SUBCELL  MODEL 

The  method  presented  here  uses  the  standard 
staggered  grid  scheme  but  now  sub-divides  or  splits  the 
existing  cell  into  8 smaller  cells  as  shown  in  Fig.  1. 
Each  field  point  may  now  be  viewed  as  a crosspoint 
surrounded  by  8 subcells  with  differing  parameters  as 
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opposed  to  being  at  the  center  of  a material  cell.  The 
material  cells  no  longer  overlap  in  space  for  the 
staggered  field  points.  A sphere,  for  example,  cutting 
through  the  standard  Yee  cell  still  is  approximated  with 
a staircase  but  with  twice  the  number  of  cells.  Unlike 
some  previous  approaches  there  are  no  special 
calculations  related  to  the  specific  geometry  modeled. 

Applying  Ampere’s  Law  around  a contour  Cx  one 
obtains: 


direction  as  shown  in  Fig.  2.  The  dielectric  sphere  has  a 
radius  of  4.5  cm.  In  this  work  we  consider  scattering  at 
two  frequencies  and  two  dielectric  constants.  Case  I is  at 

a frequency  of  2.5  GHz  with  £ — 4s0  in  the  sphere  and 

Case  II  is  at  1.768  GHz  with£  = &£0  in  the  sphere.  In 

each  case  the  wavelength  in  the  sphere  is  the  same.  The 
Mie  series  solution  for  the  fields  inside  and  outside  the 
sphere  is  compared  to  numerical  FDTD  results. 


d_ 

dt 


jD-dS^jH-d'l j. 

Si  G 


(2) 


From  Fig.  1,  assuming  that  Dz  is  constant  over  a 

surface  patch  Sx  (center  x-y  plane  in  Fig.  1)  and  is 
equal  to  the  value  at  the  center  of  the  patch,  the  surface 
integral  becomes: 


\DJS  = e O',  j,  k)Ez  | iJk  AxAy . (3) 

s, 

Assume  also  that  the  value  of  the  H field  along  each 
contour  leg  is  constant  and  equal  to  the  value  at  the 
midpoint  of  the  leg.  Since  the  material  cell  is  now  split, 
further  assume  that  at  the  cell  center  (?,  j\  k)  or 
crosspoint  of  the  8 material  cells,  that  the  effective 
dielectric  constant  is  the  average  of  the  8 surrounding 
material  cells  or : 

8 

e(i,j,k)  = e(i,j,k)  = l/8^X  . (4) 

n=\ 


Using  central-difference  approximations  for  the  time 
derivatives  (interleaved  in  time)  the  following  step- 
ahead  equation  results  - valid  for  non-cubical  cells. 


17/7+1 

ijjk 


i,jjc 


At 

+3 

£(i,j,k) 


ij-]/2Jc 
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At  y 
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u+yu 


l_ 

Ax 


i-ytjji 


(5) 


Analogous  equations  can  be  developed  for  Ey  and  Ex  . 


III.  BENCHMARK  PROBLEM  - DIELECTRIC 
SPHERE 

The  problem  considered  is  a plane  wave  scattering  from 
a dielectric  sphere.  The  incident  wave  has  the  electric 
field  polarized  in  the  z direction  and  travels  in  the  +y 


A.  Computational  Modeling 

Numerical  results  are  generated  using  a 3D  FDTD 
computational  code  that  implements  the  standard  Yee 
method  with  and  without  the  subcell  averaging.  The 
computational  grid  is  512  x 256  x 512.  At  the 
extremities  of  the  grid  Liao  second  order  absorbing 
boundary  conditions  (ABCS)  are  employed.  A sine 
wave  source  is  fed  in  at  y=5.  The  incident  wave  is 
planar  over  1/2  of  the  grid's  width  or  from  x=  128  to 
x=384,  with  a Gaussian  taper  toward  the  absorbing 
boundaries  This  finite-width  waveform  approximates 
plane  wave  conditions  in  the  vicinity  of  the  sphere.  The 


spatial  discretization  Ax 


is  set  to 


dielectric  media)  for  both  Cases  I and  II.  The  sphere  has 
a radius  of  1 5 Ax  ;however,  with  the  split  cell  model, 
the  sphere  has  a radius  of  30  subcells.  The  FDTD 
simulation  is  run  in  continuous  wave  (cw)  mode,  long 
enough  to  achieve  steady-state  conditions  over  the 
extent  of  the  grid.  The  amplitudes  of  the  vector 
components  of  the  E field  are  saved  along  chosen  line 
cuts  (along  y)  though  the  sphere.  The  program,  which 
updates  large  3D  grids,  is  written  with  OpenMP 
directives  to  run  efficiently  on  an  SGI  Origin  parallel 
computer. 


B.  Analytical  Solution 

The  Mie  solution  for  a plane  harmonic  wave 
scattering  from  a dielectric  sphere  is  well  known  and 
may  be  found  in  several  references  [6,7].  This  exact 
solution  is  a series  for  the  scattered  field  exterior  to  the 
sphere  and  the  transmitted  or  total  field  interior  to  the 
sphere.  The  incident  field  is  added  to  the  scattered  field 
to  obtain  the  total  field  outside  the  sphere.  The  resulting 
vector  components,  in  spherical  polar  coordinates,  are 
transformed  into  Cartesian  vector  components  along 
various  cuts  though  the  sphere.  This  permits  a direct 
comparison  with  the  FDTD  solution.  Care  must  be 
taken  to  insure  that  the  FDTD  cuts  correspond  exactly 
in  space  with  the  analytical  ones. 


C.  Results  and  Comparison 

Fig.  3 shows  comparisons  for  Case  I between  the 
Mie  solution  (exact),  the  8 subcell  or  split  cell  model, 
and  standard  cell  model  for  the  magnitude  of  Ez  along  a 
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straight  line  through  the  sphere.  For  cuts  parallel  to  the 
y-axis,  Ez(0,y,0)  in  Fig.  3(a)  and  Ez(8,y,0)  in  Fig.  3(b), 
the  standard  cell,  the  split  cell  and  exact  results  are  close 
and  difficult  to  distinguish.  However,  the  standard  cell 
model  exhibits  a local  amplitude  overshoot  at  the  back 
of  the  sphere.  The  subcell  model  produces  more 
accurate  field  values  at  the  back  of  the  sphere. 

Fig.  4 shows  the  comparison  between  the  Mie 
solution,  the  subcell  model,  and  the  standard  model  for 
the  magnitude  of  Ey  Along  the  cut  Ey(0,y,l)  shown  in 
Fig.4(a),  the  standard  cell  model  gives  a large 
undershoot  and  spatial  phase  error  near  the  back  of  the 
sphere,  where  there  is  a discontinuity  in  E>,  The  subcell 
model  does  not  exhibit  this  large  undershoot  or  spatial 
error.  The  standard  cell  model  also  overshoots  the  exact 
solution  near  the  front  of  the  sphere.  Fig.  4(b)  shows  the 
cut  Ey(0,y,8).  The  standard  cell  and  split  cell  are  closer 
in  amplitude  away  from  the  center  (z=0)  plane. 

Figs.  5 and  6 show  the  analogous  comparisons  for 
Case  II.  A rather  large  overshoot  is  apparent  in  Fig.  5(b) 
near  the  back  of  the  sphere  in  the  z=8  plane  for  the 
standard  cell  case.  Fig.  6(a)  best  illustrates  the  standard 
cell  errors  and  shows  large  undershoots  at  both  the  front 
and  back  of  the  sphere,  where  there  are  jump 
discontinuities  in  Ey  An  overshoot  at  the  back  of  the 
sphere  for  the  standard  cell  method  is  also  apparent  in 
Fig.  6(b).  The  split  cell  model  agrees  better  with  the 
exact  at  both  the  front  and  back  of  the  sphere. 

To  compare  the  split  cell  and  averaging  method 
with  the  similar  dielectric  subgrid  resolution  (DSR) 
method  [5],  Case  II  was  computed  using  a grid  spacing 
of  1/10  of  a wavelength  in  the  sphere.  This  coarse  grid 
accentuates  any  differences  between  the  methods.  The 
Ey  component  across  the  interface  best  illustrates  these 
differences.  Fig.  7(a)  shows  that  the  DSR  method  tends 
to  underestimate  the  jump  discontinuities  as  compared 
to  the  simpler  averaging  method.  However,  at  1/20  of  a 
wavelength  in  the  sphere,  the  DSR  and  simple 
averaging  curves  are  both  converging  to  the  Mie  result, 
as  shown  in  Fig. 7(b). 

IV.  PLANAR  INTERFACES 

The  performance  of  the  split  cell  FDTD  method  is 
further  demonstrated  by  computing  the  reflection  and 
transmission  of  a plane  wave  Gaussian  pulse  polarized 
in  the  z direction  from  a planar  interface.  The  grids  are 
chosen  to  be  128x64x128.  An  interface  tilted  at  14°  is 
simply  a wedge  angled  from  (0,64,z)  to  (128,128,z). 
The  grid  spacing  is  based  on  the  2/3  down  power  point 
for  the  pulse.  At  this  frequency  the  grid  spacing  is 
chosen  to  be  1/10  of  a wavelength  in  the  dielectric.  At 
smaller  grid  spacings  the  differences  between  the  split 
cell  and  standard  cell  become  minimal  for  a single 
pulse.  For  1/10  of  a wavelength  grid  spacing  in  the 
dielectric  a second-order  accurate  FDTD  formulation  is 
used.  The  same  cases  are  repeated  with  a coarser  grid 


spacing  of  1/7.5  of  a wavelength  but  using  a fourth- 
order  (second-order  in  time  and  fourth-order  in  space  or 
2,4  scheme)  FDTD  formulation. 

The  peak  amplitudes  of  the  transmitted  pulse  and 
reflected  pulse  are  compared  with  the  exact  amplitudes. 
In  order  to  capture  the  reflected  amplitude  the  incident 
pulse  is  subtracted  out  on  the  reflecting  side  of  the 
interface.  The  results  are  summarized  in  Tables  I-IV. 
The  percent  errors  are  shown  in  parentheses.  The  data 
show  that  the  split  cell  formulation  is  consistently  more 
accurate  for  interfaces  angled  at  14°.  The  cell  splitting 
and  averaging  is  particularly  useful  on  coarse  grids 
where  staircasing  errors  may  be  more  of  a concern.  The 
split  cell  and  standard  cell  give  the  same  transmitted  and 
reflected  amplitudes  for  0°  (normal  incidence)  thus 
illustrating  that  the  averaging  only  has  an  effect  for 
interfaces  at  angles. 

V.  CONCLUSIONS 

The  split  cell  and  averaging  method  is  shown  to 
give  improved  accuracy  across  angled  interfaces  and 
along  curved  interfaces,  such  as  a sphere.  The  method  is 
especially  useful  where  coarser  grids  are  required.  The 
technique  is  very  simple,  only  requiring  that  the 
material  cells  have  twice  the  resolution  of  the  field  cells. 
Since  the  averaging  of  the  dielectric  constants  is  done 
prior  to  the  time  stepping  of  the  fields  there  is  little 
added  computational  cost. 
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Fig.  1 Diagram  of  Yee  cell  with  8 subcells  along  with  the  position  of  the  E and  H vector  components. 
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Fig.  2 Diagram  of  incident  wave,  dielectric  sphere,  and  coordinate  system. 
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Figure  4a  and  4b  - Magnitude  of  Ey.  parallel  to  y-axis  through  (x=0,z=l)  and  (x=0,z=8).  Case  I,  f=2.5GHz  and  8=4. 
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Figures  6a  and  6b  - Magnitude  of  Ey  parallel  to  y-axis  through  (x=0,z=0)  and  (x=0,z=8).  Sphere  is  located  between 
y=l  13  and  y=143.  Case  II,  f=l  .767  GHz  and  e=8. 
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Figure  7a  - Magnitude  of  Ey  parallel  to  y-axis  through  (x=0,z=l).  Compared  are  DSR  and  simple  averaging.  Case  II, 
f=  1.767  GHz,  and  e=8,  10  points/wavelength.  Figure  8b  is  same  as  8a  but  with  20  points/wavelength. 
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Table  I 

Amplitudes  of  Transmitted  and  Reflected  Waves 
Second-Order,  £=4.0,  14°,  A/1 0 


Split  Cell 

Standard  Cell 

Exact 

Transmitted 

Amplitude 

.660 

(.3%) 

.655 

(1.06%) 

.662 

Reflected 

Amplitude 

-.335 

(3.7%) 

-.344 

(6.5%) 

-.323 

Table  II 

Amplitudes  of  Transmitted  and  Reflected  Waves 
Second-Order,  £=8.0,  14°,  A/10 


Split  Cell 

Standard  Cell 

Exact 

Transmitted 

Amplitude 

.519 

(0.0%) 

.513 

(1.16%) 

.519 

Reflected 

Amplitude 

-.477 

(1.9%) 

-.484 

(3.4%) 

-.468 

Table  III 

Amplitudes  of  Transmitted  and  Reflected  Waves 
Fourth-Order,  £=4.0,  14°,  A/7.5 


Split  Cell 

Standard  Cell 

Exact 

Transmitted 

Amplitude 

.660 

(0.3%) 

.643 

(2.8%) 

.662 

Reflected 

Amplitude 

-.325 

(.6%) 

-.351 

(8.7%) 

-.323 

Table  IV 

Amplitudes  of  Transmitted  and  Reflected  Waves 
Fourth-Order,  £=8.0,  14°,  A/7.5 


Split  Cell 

Standard  Cell 

Exact 

Transmitted 

Amplitude 

.518 

(0.3%) 

.502 

(3.3%) 

.519 

Reflected 

Amplitude 

-.468 

(0.0%) 

-.494 

(5.6%) 

-.468 

